 rm(list = ls())


 library(car)  # needed for recode function
 library(Zelig)
 library(Amelia)

#setwd("")

load("d_RUS_FIRMS_2010_MI.RData")


library(sandwich) 
library(lmtest) 


######  LAWYERS2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(lawyers2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.lawyers2<-combined.results$q.mi[2]
se3.lawyers2<-combined.results$se.mi[2]

######  COURTS2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(courts2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.courts2 <-combined.results$q.mi[2]
se3.courts2 <-combined.results$se.mi[2]

######  LAW ENF2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(lawenf2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.lawenf2 <-combined.results$q.mi[2]
se3.lawenf2 <-combined.results$se.mi[2]

######  BUR2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(bur2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.bur2 <-combined.results$q.mi[2]
se3.bur2 <-combined.results$se.mi[2]

######  COURTS_INF2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(courts_inf2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.courts_inf2 <-combined.results$q.mi[2]
se3.courts_inf2 <-combined.results$se.mi[2]

######  LAWENF_INF2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(lawenf_inf2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.lawenf_inf2 <-combined.results$q.mi[2]
se3.lawenf_inf2 <-combined.results$se.mi[2]


######  BUR_INF2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(bur_inf2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.bur_inf2 <-combined.results$q.mi[2]
se3.bur_inf2 <-combined.results$se.mi[2]



######  PSS2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(pss2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.pss2 <-combined.results$q.mi[2]
se3.pss2 <-combined.results$se.mi[2]


######  PSA2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(psa2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.psa2 <-combined.results$q.mi[2]
se3.psa2 <-combined.results$se.mi[2]

######  MAFIA2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(mafia2~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.mafia2 <-combined.results$q.mi[2]
se3.mafia2 <-combined.results$se.mi[2]

####################################


######  LAWYERS

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(lawyers~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)



coef3.lawyers<-combined.results$q.mi[2]
se3.lawyers<-combined.results$se.mi[2]

######  COURTS

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(courts~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.courts <-combined.results$q.mi[2]
se3.courts <-combined.results$se.mi[2]

######  LAW ENF2

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(lawenf~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.lawenf <-combined.results$q.mi[2]
se3.lawenf <-combined.results$se.mi[2]

######  BUR

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(bur~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.bur <-combined.results$q.mi[2]
se3.bur <-combined.results$se.mi[2]

######  COURTS_INF

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(courts_inf~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.courts_inf <-combined.results$q.mi[2]
se3.courts_inf <-combined.results$se.mi[2]

######  LAWENF_INF

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(lawenf_inf~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.lawenf_inf <-combined.results$q.mi[2]
se3.lawenf_inf <-combined.results$se.mi[2]


######  BUR_INF

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(bur_inf~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.bur_inf <-combined.results$q.mi[2]
se3.bur_inf <-combined.results$se.mi[2]



######  PSS

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(pss~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.pss <-combined.results$q.mi[2]
se3.pss <-combined.results$se.mi[2]


######  PSA

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(psa~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.psa <-combined.results$q.mi[2]
se3.psa <-combined.results$se.mi[2]

######  MAFIA

b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out<-lm(mafia~tax1_dich_r+other_firms2_r+privatized*consolidated+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data=RF$imputations[[i]])
JGM.out.rob<-coeftest(JGM.out, vcov=vcovHC(JGM.out, type="HC1"))
b.out<-rbind(b.out,JGM.out.rob[,1])
se.out<-rbind(se.out,JGM.out.rob[,2])
}
combined.results <- mi.meld(q = b.out, se = se.out)

coef3.mafia <-combined.results$q.mi[2]
se3.mafia <-combined.results$se.mi[2]


#########
#FIGURE 2 CODE
#########

# will produce PDF output

#pdf("fig_2", height =6, width = 9) #open pdf device



layout(matrix(c(3,1,2), 1,3),widths=c(1,3,2.4))  ## allows for three graphs, in matrix function is "plot1" through "plot 3"; then commands for 1 rows, 3 columns; c(3,1,2) puts 3rd plot first; widths sets relative width of columns

#PLOT 1

#Create Vectors for coefs and standard errors for each model, and variable names
   
coef.vec.1 <-c(coef3.mafia2,coef3.psa2,coef3.pss2,coef3.courts_inf2,coef3.lawenf_inf2,coef3.bur_inf2,coef3.lawyers2,coef3.courts2  ,coef3.lawenf2,coef3.bur2)
se.vec.1 <- c(se3.mafia2,se3.psa2,se3.pss2,se3.courts_inf2,se3.lawenf_inf2,se3.bur_inf2,se3.lawyers2,se3.courts2,se3.lawenf2,se3.bur2)

 
var.names <- c("criminal racket", "private security \nagency ","internal security \nservice ", "courts \n(with connections)","law enforcement \n(unofficial)","bureaucrats \n(unofficial)","lawyers \n(out of court)","courts", "law enforcement \n(official)", "bureaucrats \n(official)")

y.axis <- length(var.names):1#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis
adjust <- .18 #create object that we will use to adjust points and lines up and down to distinguish between models
    

par(mar=c(3,8,7,1), lheight = .8)#set margins for regression plot, c(bottom, left, top, right)
plot(coef.vec.1, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.3, #plot model 1 coefs using black points (pch = 19, default = black), adding the "adjust amount" to the y.axis indicator to move points up
 
 
xlim = c(-2.1, 2.1), xaxs = "r", #set limits of x-axis so that they include mins and maxs of        
ylim = c(min(y.axis)-.2, max(y.axis)+.2),yaxs="r", main = "",cex.main=1)
   title("Property Dispute",line=5.8,cex.main=1.7)
axis(1,at = seq(-4,4, by = 1), label = seq(-4,4, by = 1), mgp = c(0,1,0), cex.axis = 1.2)#add x-axis and labels; c(title,labels,line) "pretty" creates a sequence of  equally spaced nice values that cover the range of the values in 'x'-- in this case, integers
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.2)#add y-axis and labels; las = 1 makes labels perpendicular to y-axis
#axis(3,pretty(coef.vec.1, 3))#same as x-axis, but on top axis 
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")#draw light dotted line at each variable for dotplot effect
box(bty="l")#draw box around plot
segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis, coef.vec.1+qnorm(.975)*se.vec.1, y.axis, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.1-qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1-qnorm(.95)*se.vec.1, y.axis +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.1+qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1+qnorm(.95)*se.vec.1, y.axis+.035, lwd = 1.1)   #confidence intervals
abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing
points(coef.vec.1, y.axis, pch = 21, cex = 1.3, bg = "white" ) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color


arrows(.2,10.8,1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)
arrows(-.2,10.8,-1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)

text(.2,11.45, "tax violators more \nlikely to use \nthan tax compliers",xpd=TRUE,adj=0,cex=1.2)
text(-.2,11.45, "tax violators less \nlikely to use \nthan tax compliers",xpd=TRUE,adj=1,cex=1.2)

                
#### PLOT 2

coef.vec.1 <- c(coef3.mafia,coef3.psa,coef3.pss,coef3.courts_inf,coef3.lawenf_inf,coef3.bur_inf,coef3.lawyers,coef3.courts  ,coef3.lawenf,coef3.bur)
 se.vec.1 <- c(se3.mafia,se3.psa,se3.pss,se3.courts_inf,se3.lawenf_inf,se3.bur_inf,se3.lawyers,se3.courts,se3.lawenf,se3.bur)

y.axis <- length(var.names):1#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis
adjust <- .18 #create object that we will use to adjust points and lines up and down to distinguish between models
    
par(mar=c(3,2,7,2), lheight = .8)#set margins for regression plot, c(bottom, left, top, right)
plot(coef.vec.1, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.3, #plot model 1 coefs using black points (pch = 19, default = black), adding the "adjust amount" to the y.axis indicator to move points up
 
 xlim = c(-2.1, 2.1), xaxs = "r", #set limits of x-axis so that they include mins and maxs of 
        #coefficients + .95% confidence intervals and plot is symmetric
 ylim = c(min(y.axis)-.2, max(y.axis)+.2),yaxs="r", main = "",cex.main=1)
   title("Debt Dispute",line=5.8,cex.main=1.7)
axis(1,at = seq(-4,4, by = 1), label = seq(-4,4, by = 1), mgp = c(0,1,0), cex.axis = 1.2)#add x-axis and labels; c(title,labels,line) "pretty" creates a sequence of  equally spaced nice values that cover the range of the values in 'x'-- in this case, integers
#axis(2, at = y.axis, label = "", las = 1, tick = T, cex.axis =1)#add y-axis and labels; las = 1 makes labels perpendicular to y-axis
#axis(3,pretty(coef.vec.1, 3))#same as x-axis, but on top axis 
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")#draw light dotted line at each variable for dotplot effect
#box(bty="l")#draw box around plot
segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis, coef.vec.1+qnorm(.975)*se.vec.1, y.axis, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.1-qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1-qnorm(.95)*se.vec.1, y.axis +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.1+qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1+qnorm(.95)*se.vec.1, y.axis +.035, lwd = 1.1)   #confidence intervals
abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing

points(coef.vec.1, y.axis, pch = 21, cex = 1.3, bg = "white" ) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color

arrows(.2,10.8,1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)
arrows(-.2,10.8,-1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)

text(.2,11.45, "tax violators more \nlikely to use \nthan tax compliers",xpd=TRUE,adj=0,cex=1.2)
text(-.2,11.45, "tax violators less \nlikely to use \nthan tax compliers",xpd=TRUE,adj=1,cex=1.2)

##########
library(pBrackets)

par(mar=c(3,.5,7,1.5))#set margins for regression plot, c(bottom, left, top, right)

#left.side <- .45#use this to manipulate how far segments are from y-axis
    #note:  getting braces and text in proper place requires much trial and error

plot(seq(0,1,length=length(var.names)), y.axis, type = "n", axes = F,  xlab = "", ylab = "")#call empty plot using type="n"

brackets(1, 1.2, 1, 4, h = NULL, ticks = 0.5, curvature = 0.5, type = 1, col = 1, lwd = 1, lty = 1, xpd = TRUE)
text(.55, 1.6, labels=expression(paste("Legal Strategies")), adj=c(0,0),xpd=TRUE,srt=90, cex=1.2)

brackets(1, 5, 1, 7, h = NULL, ticks = 0.5, curvature = 0.5, type = 1, col = 1, lwd = 1, lty = 1, xpd = TRUE)
text(.55, 5.05, labels=expression(paste("Illegal Strategies")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)
text(.72, 5.3, labels=expression(paste("(corruption)")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)

brackets(1, 7.9, 1, 9.8, h = NULL, ticks = 0.5, curvature = 0.5, type = 1, col = 1, lwd = 1, lty = 1, xpd = TRUE)
text(.55, 7.9, labels=expression(paste("Illegal Strategies")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)
text(.72, 8.3, labels=expression(paste("(violence)")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)


#dev.off()




